function [bellows_volt, flow, signal]=gernerate_representative_breath_toolbox(patient)
    display('Generating a representative breath, from patient.bellows_voltage_drifted ...')
    load([patient.folder_static '/bellows_voltage_scans_drifted.mat'])
    x =  bellows_voltage_scans_drifted;
    signal=detrend(x);
    mean_dv_db=2;%mean(data);
    
    [maxtab,mintab]=peakdet(signal(:,1),0.1*mean_dv_db);
    l=length(mintab(:,2));
    n=0;
    k=[]; %#ok<NASGU>
    for i=1:l-1
        k=signal(mintab(i,1):mintab(i+1,1),1);
        lk(i,1)=length(k); %#ok<AGROW>
        if size(k,1)>n
            n=size(k,1);
        end
    end
    %%
    %%
    % PHASE SORTING AND LINEAR DRIFT CORRECTION
    d=zeros(n,l)/0;
    for i=1:l-1
        k=signal(mintab(i,1):mintab(i+1,1),1);
        sk(i)=length(k); %#ok<AGROW,NASGU>
        hold on
        start_t=1;
        end_t=length(k);
        start_a=k(1);
        end_a=k(length(k));
        mx=polyfit([start_t,end_t],[start_a,end_a],1);
        C=polyval(mx,start_t:end_t);
        k2=k-C';
        d(1:size(k2,1),i)=k2(:);
        [maxd,mind]=peakdet(d(:,i),0.1*mean_dv_db); %#ok<NASGU>
        Td(i)=maxd(:,1); %#ok<AGROW>
        Ad(i)=maxd(:,2); %#ok<AGROW>
    end
    %%
    %%
    % Amplitude BASED GATING
    mg=zeros(length(Ad),1)/0;
    mgg=zeros(length(Ad),1)/0;
    A_gate=zeros(length(Ad),1);
    ub=mean(Ad)+std(Ad);
    lb=mean(Ad)-std(Ad);
    for jj=1:length(Ad)
        mg(jj)=Ad(jj)<ub;
        if mg(jj)<1
            A_gate(jj)=NaN;
        else
            if mg(jj)>0
                A_gate(jj)=Ad(jj);
            end
        end
    end
    for jjj=1:length(Ad)
        mgg(jjj)=Ad(jjj)>lb;
        if mgg(jjj)<1
            A_gate(jjj)=NaN;
        else
            if mgg(jjj)>0
                A_gate(jjj)=Ad(jjj);
            end
        end
    end
    for j=1:l-1
        if mg(j)==mgg(j)
            Da(:,j)=d(:,j); %#ok<AGROW>
        else
            Da(:,j)=(1:length(d(:,j)))*NaN; %#ok<AGROW>
        end
    end
    %%
    %%
    % PHASE BASED GATING
    ng=zeros(length(lk),1)/0;
    ngg=zeros(length(lk),1)/0;
    ub1=mean(lk)+std(lk);
    lb1=mean(lk)-std(lk);
    for jj=1:length(lk)
        ng(jj)=lk(jj)<ub1;
        if ng(jj)<1
            T_gate(jj)=NaN; %#ok<AGROW>
        else
            if ng(jj)>0
                T_gate(jj)=lk(jj); %#ok<AGROW>
            end
        end
    end
    for jjj=1:length(lk)
        ngg(jjj)=lk(jjj)>lb1;
        if ngg(jjj)<1
            T_gate(jjj)=NaN; %#ok<AGROW>
        else
            if ngg(jjj)>0
                T_gate(jjj)=lk(jjj); %#ok<AGROW>
            end
        end
    end
    for j=1:l-1
        if ng(j)==ngg(j)
            Dp(:,j)=d(:,j); %#ok<AGROW>
        else
            Dp(:,j)=(1:length(d(:,j)))*NaN; %#ok<AGROW>
        end
    end
    %%
    %%
    % COMBINED AMPLITUDE AND PHASE GATING
    for h=1:length(Da(1,:))
        for hh=1:length(Da(:,h))
            if Dp(hh,h)==Da(hh,h)
                Dap(hh,h)=Dp(hh,h); %#ok<AGROW>
            else
                Dap(hh,h)=NaN; %#ok<AGROW>
            end
        end
    end
    %%
    %%
    % PHASE SHIFT ALIGNMENT
    Tc=Td-round(mean(Td));
    D=padarray(Dap,max(abs(Tc)));
    D(D==0) = NaN;
    DD=zeros(length(D(:,1)),l);
    for ii=1:l-1
        DD(:,ii)=circshift(D(:,ii),-1*Tc(ii));
        ind=find(DD==0);
        DD(ind)=NaN; %#ok<FNDSB>
    end
    %%
    %%
    % NORMALIZE THE AMPLITUDE TO MAXIMUM VALUE
    DAnorm=zeros(length(D(:,1)),l);
    max_D=zeros(l,1);
    for pp=1:l
        temp1=find((~isnan(DAnorm(:,pp)))); %#ok<NASGU>
        max_D(pp)=max(DD(:,pp));
        DAnorm(:,pp)=(DD(:,pp))/(max_D(pp));%***
    end
    %%
    %%
    % FIND MEAN POSITION FOR EACH TIME POSITION
    mDAnorm=(zeros(length(DAnorm(:,1)),1))/0;
    for ppp=1:length(DAnorm(:,1))
        rel=find((~isnan(DAnorm(ppp,:))));
        temp2=DAnorm(ppp,rel); %#ok<FNDSB>
        if(isempty(temp2))
            temp_D(:,ppp) = NaN; %#ok<AGROW,NASGU>
        else
            mDAnorm(ppp,1)=mean(temp2);
        end
    end
    %%
    %%
    %INHALE
    %will be averaging along the x-axis so vect_x is now
    %the vector of output values (from the interpolation) that
    %we will be averaging over
    a=find(mDAnorm==max(mDAnorm));%==1;***
    lh=DAnorm(1:a,:);
    vect_xl=NaN(size(lh,1),size(lh,2),'double');
    vect_yl=single(0:1/(size(lh,1)-1):1);  %these are the points at which we will be evaluating our functions
    for ii=1:l
        holder=lh(find(~isnan(lh(:,ii))),ii); %#ok<FNDSB>
        offset=find(~isnan(lh(:,ii)),1,'first')-1;
        holder=holder+ 1E-5*rand(length(holder),1);
        if(~isempty(holder))
            vect_xl(:,ii)=interp1((squeeze(holder)),[1:size(holder,1)]',vect_yl','linear')+offset; %#ok<NBRAK>%%'linear'
        end
    end
    val_meanl=zeros(size(lh,1),1,'double');  %the mean of the breaths at a given y
    for ii=1:size(lh,1)
        val_meanl(ii)=mean(vect_xl(ii,find(~isnan(vect_xl(ii,:))))); %#ok<FNDSB>
    end
    clear holder offset
    %%
    %%
    %EXHALE
    %will be averaging along the x-axis so vect_x is now
    %the vector of output values (from the interpolation) that
    %we will be averaging over
    rh=DAnorm(a:length(DAnorm),:);
    vect_xr=NaN(size(rh,1),size(rh,2),'double');
    vect_yr=single(1:-1/(size(rh,1)-1):0);  %these are the points at which we will be evaluating our functions
    for ii=1:l
        holder=rh(find(~isnan(rh(:,ii))),ii); %#ok<FNDSB>
        offset= find(~isnan(rh(:,ii)),1,'first')-1;
        holder=holder+ 1E-5*rand(length(holder),1);
        if(~isempty(holder))
            vect_xr(:,ii)=interp1(squeeze(holder),[1:size(holder,1)]',vect_yr','linear') +offset; %#ok<NBRAK>
        end
    end
    val_meanr=zeros(size(rh,1),1,'double');  %the mean of the breaths at a given y
    for ii=1:size(rh,1)
        val_meanr(ii)=mean(vect_xr(ii,find(~isnan(vect_xr(ii,:))))); %#ok<FNDSB>
    end
    clear holder offset
    %%
    %%
    % PLOT FLOW RATE
    real_D=max_D(find(~isnan(max_D)));
    b=cat(1,vect_yl',vect_yr(2:end)');
    bb=cat(1,val_meanl,val_meanr(2:end)+val_meanl(end));
    bb_isnan=find(isnan(bb));
    bb(bb_isnan)=[];b(bb_isnan)=[];
    bothx=bb;
    bothy=b;
    nbothx=bothx(1):(bothx(end)-bothx(1))/(size(bothx,1)*10):bothx(end);
    nbothy=interp1(bothx,bothy,nbothx,'linear');
    clear bothx bothy
    o=diff(nbothx(1:2));
    nbothx=cat(2,nbothx(1100:end),nbothx(2:4000)+nbothx(end)-nbothx(2)+o);
    nbothy=cat(2,nbothy(1100:end),nbothy(2:4000));
    bothx=nbothx;
    bothy=nbothy;
    clear nbothy
    w=ceil(0.05*length(bothy));
    warning off %#ok<WNOFF>
    for i=w+1:(length(bothy)-w)
        index_start=max(1,i-w);
        index_end=min(length(bothy),i+w);
        index=w+min(0,i-w);
        p=polyfit((index_start:index_end),bothy(index_start:index_end),2);
        nbothy(i-w)=polyval(p,index_start+index);
    end
    clear bothy
    bothy=nbothy;
    W=ceil(0.01*length(bothy));
    for iii=2:(length(bothy)-W)
        v=[iii-1:iii+W]; %#ok<NBRAK>
        fl(iii-1)=(diff(bothy(v))/diff(bothx(v+w)))*mean(real_D); %#ok<AGROW>
    end
    flows=(zeros(length(DAnorm(:,1))-1,l))/0;
    for e=1:length(DAnorm(1,:))
        flows(:,e)=diff(DAnorm(:,e))*max_D(e);
    end
    figure;
    % subplot(2,2,1)
    hold on
    for ee=1:length(DAnorm(1,:))
        plot(DAnorm(2:length(DAnorm(:,ee)),ee)*max_D(ee),flows(:,ee)*100);
    end
    plot(bothy(2:end-W)*mean(real_D),fl*100,'r','linewidth',3)
    title(sprintf('Patient %d',patient.pt),'fontsize',14)
    xlabel('Tidal Volume [L]','fontsize',14)
    ylabel('Flow [L/s]','fontsize',14)
    set(gca,'FontSize',14)
    box off
    hold off
    xlim([0 max(signal)-min(signal)])
    bellows_volt_tmp=bothy(2:end-W)*mean(real_D);
    bellows_volt = bellows_volt_tmp + mean(x);
    flow=fl*100;
    
